Coupling vibration analysis of heat exchanger tube bundles under different stiffness conditions

A two-dimensional tube bundles fluid–structure coupling model was developed using the CFD approach, with a rigid body motion equation and the Newmark integral method. The numerical simulations were performed to determine the vibration coupling properties between various tube bundles of stiffness. Take the corner square tube bundles with a pitch ratio of 1.28 as the research object. The influence of adjacent tubes with different stiffness on the vibration of the central target tube was analyzed. The research results show that the vibration characteristic of tube bundles is affected by the flow field dominant frequency and the inherent frequency of tube bundles. The vibration of adjacent tube bundles significantly impacts the amplitude and frequency of the central target tube. The equal stiffness and large stiffness tubes upstream or downstream inhibit the vibration displacement of the target tube to some extent. The low-stiffness tubes upstream or downstream significantly enhanced the amplitude of the target tube. The findings can be used to provide a basis for reasonable design and vibration suppression of shell-and-tube heat exchangers.

www.nature.com/scientificreports/critical flow velocity of square-arranged bundles, and the calculated results agree well with the experimental data.Tang et al. 13 proposed a high-precision, low-cost CFD/CSD high-order coupling algorithm.Fluent software was combined with the dynamic mesh method to systematically study the dynamic response of the structure of the tube bundles in the dense heat exchange tube bundles.Tube drag, lift, frequency, and amplitude were analyzed in the time and frequency domain.
In a typical shell-and-tube heat exchanger, the stiffness of tubes is different under different constraints.The tubes can be generally divided into A, B, and C types according to the different number of baffle plates.As shown in Fig. 2, the failure shell and tube heat exchanger have five baffle plates.A heat exchange tube passes through three baffle plates and two baffle gap areas while the stiffness is centered.B type of heat exchange tube through all five baffle plates, the large stiffness.The C type of heat exchange tubes passes through two baffle plates and three baffle notched areas with small stiffness.
The influences of coupling vibration of tube bundles on the central target tube are investigated when the stiffnesses of adjacent tubes are different, especially for the multiple flexibly-mounted tube bundles 16,17 .As shown in Fig. 1, the location of failure heat exchange tubes is the interface area of the two types of tubes, illustrating that the stiffnesses of tube bundles in the transverse flow are quite different.It should be considered when the coupling vibration characteristics between tubes are analyzed.
Fluent software is adopted to calculate the flow between tube bundles numerically.According to the actual cases of tube bundle damage, the vibration behavior of adjacent multiple bundles and the interaction mechanism between tubes were studied considering the non-uniform bundle stiffness combination.The coupling vibration of the central target tube and adjacent tubes with different stiffness was analyzed from the aspects of the flow field, structural response, and the movement path of the tube bundles to reveal the coupling vibration mechanism between non-uniform bundle stiffness.

Turbulence model
Due to the small gaps between tubes, the fluid velocity in the main flow area significantly differs from the velocity between the gaps.Meanwhile, the flow around produces a large number of vortex structures.Therefore, it is essential to select a suitable turbulence model to simulate the flowfield between dense tube bundles to describe the interaction between the fluid and structure honestly.
Researchers have used different turbulence models to numerically calculate the flow field induced by tube beam vibration.For example, Bao 19 used RNG k-ε to calculate the three-dimensional flow field between tube bundles with a pitch ratio 1.28.By comparing with the experimental data, Zhao 20 believed that the Transition SST model is more suitable for numerical simulation of turbulent flow field between close-packed tube arrays.Darwish 21 stated that the SST k-ω was more reliable in solving the flow between a rotated square tube array.
The turbulence model selected will be verified by comparing the computational results of different turbulence models with experimental data.Based on the velocity field experimental data of the staggered bundles of Balabani 22 , as shown in Fig. 3, a two-dimensional calculation model was established.The calculation mesh is shown in Fig. 4.There are 108 nodes on the tube bundle boundary.The number of encryption layers was 6, the (1)  www.nature.com/scientificreports/initial height was 0.2 mm, and the mesh growth rate was 1.1.The RNG k-ε, SST k-ω, and Transition SST model were adapted for numerical simulation, and the velocity field distribution was obtained.The comparison of x velocities calculated by different turbulence models is shown in Fig. 5.Total of 4 positions, x/d = 0, 2.5, 4.2, and 6.7, are compared, respectively.Maximum relative error of different turbulence models are shown in Table 1.At the initial position, all three models agree well with the measured results.The transition SST has the most significant calculation error, but the maximum relative error is only 6.3%.At x/d = 2.5, the overall trend of the calculated results of the three models is consistent with the experimental result.Still, in numerical terms, the SST k-ω model is the most compatible with the experiment results, with a maximum relative error of 13.8%.At the same time, the RNG k-ε and the transition SST have the maximum relative error of 27.6% and 38.7%, respectively.With flow development, the difference between different models is more significant.At x/d = 4.2, the maximum relative error of transition SST is 36.3%.The SST k-ω is also the closest to the actual measurement value, with a maximum error of 16.2%.At x/d = 6.7, the maximum error of SST k-ω is 9.3%.According to the above results, the error of the SST k-ω model is the smallest, and it is in the best agreement with the measured results.Therefore, the SST k-ω model is adopted in this work.
The SST k-ω model was proposed by Menter 23 .Compared to most RANS models, the SST k-ω model provides better prediction of flow separation and performs well at lower pressure gradients, with high accuracy and economy.Different from standard k-ω, the SST k-ω model modifies the model constants and takes into account the transfer of turbulent shear stress, which uses k-ω in the inner region of the boundary layer and switches to k-ε in the free shear flow 24 .where P = τ ij ∂u j ∂x j , k is turbulence kinetic energy.ω is turbulence-specific dissipation rate.µ t is turbulence viscosity coefficient.

Dynamics modeling of tube bundles vibration
The fluid-induced vibration of the tube bundles involves the interaction between the fluid and solid domain, which is challenging to calculate numerically 25 .Based on the computational fluid dynamics method, this paper establishes a simplified fluid-tube coupling vibration calculation model using the dynamic mesh method.

Simplification of vibration tube bundle structure model
Because the calculation model that directly considers the coupling of tube bundle deformation and flow field change requires a lot of calculation and many iterations, the accuracy is difficult to guarantee.As shown in Fig. 6, the vibration model of the tube bundle is simplified into a two-dimensional spring-damp-mass system to improve the calculation efficiency 26,27 .
The stiffness, damping, and mass of the simplified tube bundles were calculated by referring to references 15,28,29 , respectively, and the values are shown in Table 3.

Solution of rigid body motion of bundle
The fluid force on the tube bundle is obtained by solving the fluid motion equation.Taking the fluid force as known, the rigid body motion equation on the tube bundle can be solved, and the fluid-structure coupling of the tube bundle system is realized.The rigid body's motion equation at the bundle's boundary is solved by the Newmark-β integral method 30,31 .
The vibration differential equation of the tube is where m, c, k, and F are mass, damping, stiffness of the spring-damp-mass system, and fluid force applied to the tube bundles.x , ẋ , ẍ are the displacement, velocity and acceleration of tube 32,33 .
According to the Newmark-β method, where x t , ẋt and ẍt are the displacement, velocity, and acceleration at t. x t+ t , ẋt+ t and ẍt+ t are the displace- ment, velocity and acceleration at t + Δt, γ and β are coefficient of Newmark-β method.Δt is the time step.
(3)  9) and ( 10), Then, the vibration differential equation of the system at time t + Δt is x t+ t can be expressed as K * x t+ t = F * t+ t .Where where K * is effective stiffness.F * t+ t is effective load.When the mass, stiffness, and damping of the vibrating tube bundle are obtained, the dynamic response of the next time can be obtained by substituting the previous time's dynamic response into the system's initial state.

Calculation domain and model
Due to the large number of tubes in the heat exchanger, the cost of numerical calculation is significant, making it unrealistic to model and study the vibration of all tubes.It makes sense to reduce the number of tubes through reasonable simplification 34,35 .The three-dimensional model of the tube bundle is simplified into a  The geometric structure of the numerical calculation model is shown in Fig. 7. Take the corner square arrangement tube bundle as the research object.The heat exchange tube is simplified into a bundle with seven rows and columns.Fluid flows in from the left side and out from the right, and the top and bottom sides are symmetric.The distance from the inlet to the first tube bundle is ten times the tube bundle diameter, and from the last tube bundle to the outlet is 20 times the tube diameter.The tube bundle has a diameter of 25 mm and a pitch ratio of 1.28.
Figure 7 also gives the scheme of the tube bundle stiffness combination.Tube 0 is the target tube, and tubes 1-8 are flexible tubes, which have been marked in the dotted box.The calculation method is the same as the target tube, and self-excited vibration can occur.The unlabeled tubes outside the dotted box are rigid, assuming no vibration occurs in the flow field and no particular definition is required.Three kinds of stiffness were considered: large, normal, and small.Large stiffness was taken as two times the normal stiffness, and the small stiffness was 1/2 of the normal stiffness.In the calculation, only the influence of stiffness change is considered.Thus, the mass of the tube bundle remains unchanged, and the inherent frequency changes due to the stiffness change.The stiffness and corresponding inherent frequency are shown in Table 4.
The stiffness of tube 0 is kept unchanged, and the tube bundles' stiffness is changed upstream or downstream of the target tube.The actual structure of the heat exchanger is referred to ensure that there are only two stiffnesses in the flexible tube bundle.The upstream and downstream are not repeated to determine the stiffness combination scheme.
Computational mesh and mesh independence study Mesh independence study is very important for the numerical study authenticity.Based on the number of mesh cells, different mesh densities are used to find a compromise between computational cost and accuracy.Numerical calculations were performed on meshes of 140,000, 180,000, 250,000, 310,000 and 400,000 (corresponding to mesh spacings of 1 mm, 0.8 mm, 0.6 mm, 0.5 mm and 0.4 mm) with inlet velocity of 0.5 m/s.To simplify the calculations, all tube bundles are set to be rigid.The lift coefficients of the target tube (tube 0) are compared under different meshes to determine the appropriate mesh size.The results are shown in Fig. 8.After the number of meshes reaches 310,000, as shown in Fig. 8, the calculation results are almost unaffected by the sparsity of the grids.Relative error between number of 310,000 and 400,000 is just 2.7%.So this paper's research is conducted using 310,000 grid cells.
The computational mesh is shown in Fig. 9. Due to the calculation requirements of the dynamic mesh, the combining structure and unstructured mesh is used, in which an unstructured mesh was used in the tube bundles gaps.The near tube wall zones were divided into a structure mesh.The wall surface of the bundle was encrypted.The number of encryption layers was 8, the initial height was 0.2 mm, and the mesh growth rate was 1.1.The bundle's front and rear flow areas adopt a quadrilateral mesh structure.

Boundary condition and solution settings
Because of the need to simultaneously calculate the vibration characteristics of the tube bundle and realize coupling with the fluid, it is necessary to perform transient calculation 36 .The boundary condition and solution Settings are shown in Table 5.

Experimental facility
A vibration experimental facility for a single flexible tube bundle is established to validate the fluid-structure interaction computational model based on dynamic mesh and linear solution methods.The experimental facility is shown in Fig. 10.It consists of three parts: a water circulation system, a data acquisition and analysis system, and the tube bundle vibration unit.The water circulation system comprises a pump, a supply adjustment system, flow and pressure measuring devices, a water tank, and pipe.Within the tube bundle vibration testing unit, the internal space measures 130 mm × 130 mm in cross-section.The experimental tube bundle is arranged in a 3 × 3 square configuration, with a flexible tube at the center, while the other tubes are firmly connected to the tube plate.Springs are employed to constrain the middle tube bundle to the fixed tube plate.Eight springs are evenly distributed at four positions-horizontally and vertically-at both ends of the flexible tube.Sensors can   Time step 0.0001 s be fixed and connected through bolts at the ends of the central flexible tube.The data acquisition and analysis system is the DH5981 data acquisition system, using the complementary DHDAS software for signal analysis and processing.The sensors utilized are IEPE underwater piezoelectric acceleration sensors, with their relevant parameters detailed in Table 6.

The comparison between numerical computational results and experimental data
In Figs.11 and 12, the acceleration time course curves and spectral diagrams for both experimental and numerical simulations are presented under the condition of an inlet flow velocity of 0.1 m/s.It can be observed that, overall, the trends of the two are quite similar, with the simulated values slightly lower in magnitude compared to the experimental values.In the frequency spectrum, the main frequencies of the two are also relatively close.

Measurement systems
Data acquisition instrument Signal analysis and processing system Data acquisition and analysis system Water circulation system

Analysis of coupling vibration characteristics of tube bundles Vorticity distribution between tube bundles
Figure 13 shows the vorticity distribution between multiple flexibly-mounted tubes with different stiffness.SFM represents a single flexible mounted tube bundle with A target tube stiffness.The meanings of other labels are shown in Fig. 7.The vorticity distribution shows the process of vortex formation, shedding, and adhesion between the tube bundles 37 .The first tube bundles generate the vortices, which shed, attach, and move downstream along the flow direction.Due to the interaction between the tightly packed tube bundles, the vortices between the bundles could not fully develop until the last row of tube bundles.The vortices with different shedding characteristics fell off and formed a complete Karman vortex street 38 .
When the flexible tubes are added around the target tube, the vortex shape remains unchanged, the vortex separation point is more significant than 120°, and the vortex width is narrow.In general, the overall vortex strength of the multi-flexible tube bundles is lower than the single flexibly-mounted tube bundles, so the pressure fluctuation caused by the vortex shedding will also decrease.The vortex shedding patterns on the nine flexible mounted tubes are basically the same for AAA, BAA, and AAB.The influence of the upstream tube on vortex generation and shedding of the downstream tube is weakened, and the parallel vortex on the target tube is less than that on a single flexible mounted tube.When the upstream or downstream are large stiffness tubes, the vibration law of the target tube tends to a single flexible mounted tube because the vibration amplitude of a large stiffness tube decreases.
As a comparison, the upstream or dowstream are small stiffness for CAA and AAC.Under the same condition, the amplitude of the small stiffness tube is larger.The displacement of tube 1 in the lift direction is large, resulting in a certain distance between tube 1 and downstream tube 0 in the lift direction.When the red vortex at the lower part of tube 1 falls off and is attached to tube 0 again, it is already below the former freeze point of tube 0, and the pressure at the lower part of tube 0 increases.Meanwhile, the blue vortex at the upper part of tube 0 is in the shedding stage, and the upper pressure further decreases.
Combining the two causes a larger amplitude of the target tube 0. For AAC, when the low stiffness tube is added downstream, the large vibration of 4 or 6 tubes prevents the shedding of the vortex of 0 tubes to some extent, the pressure release is blocked, and the surface force of 0 tubes is increased.When the vortices above

Time-frequency analysis of tube bundle vibration
The amplitude of target tube 0 is analyzed by time-frequency.The result is shown in Fig. 14.From the amplitude time history curve, and it can be found that compared with the single flexible mounted tube, the amplitude of the target tube is smaller under the equal stiffness AAA, the upstream or downstream large stiffness BAA or AAB.The amplitude of the target tube is significantly suppressed.It should be noted that the upstream and downstream tube bundles are moderately stimulated under this working condition, and the vibration amplitude and phase of the upstream and downstream tubes are in a suitable range, which produces a good vibration suppression effect.For CAA and AAC, the amplitude of the target tube is equal to or larger than that of the single flexible tube, which can be considered to be caused by the flow block and vortex detachment deviation caused by the large motion of the low-stiffness tube.In conjunction with Table 7, it can be seen that the amplitude of AAA is the smallest, and its root mean square amplitude is reduced by more than 50% compared with a single flexible tube.The amplitude of BAA and AAB is about 10% higher than AAA's.The amplitude of the target tube is larger for CAA and AAC, with small stiffness.
From the spectrum diagram, we can find that the dominant frequency of the vibration of the target tube is about 25, which is the vortex-shedding frequency of the flow field.The addition of upstream and downstream flexible tubes complicates the target tube's amplitude frequency 39 .In addition to the primary frequency, there are multiple peaks, which should include the inherent frequency of the bundle and the exciting frequency of turbulence.

Trajectories of each flexible mounted tube under different stiffness combinations
The motion trajectories of all flexible tubes were analyzed, and the motion trajectories of tube bundles within 0.3 s were selected, as shown in Figs. 15, 16, 17, 18, 19 and 20.For the equal stiffness AAA, each tube presents a typical eight or elliptic vibration trajectory 40 , and the amplitude of the tube bundle in the lift direction is slightly larger than that in the drag direction.The tube bundles move approximately periodically.Among them, the amplitudes of the upstream tube bundles 1, 2, 3, 7, and 8 are slightly larger than that of the downstream tube bundles 4, 5, and 6.The front of the upstream tube bundles are direct flow or fixed tubes, and the coupling characteristics        www.nature.com/scientificreports/ between other tubes and a single flexible tube have similarities, so the amplitude is more prominent.The front flexible tubes affect the downstream tube bundles, and the vibration is suppressed.
For BAA, as shown in Fig. 17, the upstream tubes 1, 2, and 8 are large stiffness tubes.Under the premise of no resonance, the amplitude of large stiffness tubes is significantly lower than that of small stiffness tubes, as reflected in the trajectory.The amplitude of large stiffness tubes is one order of magnitude smaller than that of normal stiffness tubes.According to the trajectory, the amplitude of the target tube and other normal stiffness tube bundles downstream did not decrease after adding the upstream large stiffness tube, and the amplitude of the lift direction was more significant than the resistance direction.
The downstream tubes 4, 5, and 6 for AAB are large stiffness tubes.As shown in Fig. 18, the amplitude of the large stiffness tube is significantly lower than that of the normal stiffness tube, regardless of whether it is upstream or downstream.Adding the downstream tubes with large stiffness increases the amplitude of the target tube slightly.The influence of stiffness change of the downstream tube bundle on the target tube is more potent than that of the upstream tube bundle.
As shown in Fig. 19, when small stiffness tubes are added upstream, all flexible tube bundles amplify.Among them, the displacement of small stiffness tubes 1, 2, and 8 in the lift direction is about an order of magnitude larger than that of normal stiffness tubes, and the maximum amplitude on one side is greater than 2 mm, that is, more than 0.08 days.The three tubes with this stiffness are already in resonance from the trajectory.Because of the fluid excitation force generated by the larger amplitude of the upstream small stiffness tubes, the vibration of the downstream normal stiffness tubes is also enhanced, and the increase of the 3 and 7 tubes is the largest.
The AAC is a downstream small stiffness tube bundle.As shown in Fig. 20, under this condition, the small stiffness tubes 4, 5, and 6 have large amplitudes, and the maximum amplitude on one side is more significant than 2.3 mm, which is about 0.1 days.These tubes should be in a resonance state.The trajectories of the three small-stiffness tubes are also slightly different.Tube 5 shows noticeable vibration enhancement in the lift direction.However, the resistance direction is the same as that of the normal stiffness tubes, and the obstruction of the Tube 4 and the Tube 6 in the flow direction.Because of the influence of the small stiffness tube downstream, the vibration of the normal stiffness tube bundle upstream is enhanced.

Interactions between pipes at different flow rates.
The vibration response of the AAA combination is analyzed at different flow velocities, and the amplitude-time curves and frequency-domain curves are presented in Fig. 21 The amplitude and dominant frequency of the vibration of target tube 0 both increase with the increase in flow velocity.Additionally, non-dominant frequency components in the spectrum also increase, including both the inherent frequencies of the tube bundle and indistinct turbulent disturbance components.This suggests that the increase in flow velocity enhances the coupling vibration effects between the tube bundles and the degree of turbulence-induced vibration.
Figure 22 illustrates the maximum amplitude and root-mean-square amplitude of target tube 0 concerning various stiffness combinations at different incoming flow velocities.Under AAA, AAB or BAA, the target tube's amplitude notably decreases compared to configurations with lower stiffness in the upstream or downstream sections.Particularly, the upstream higher stiffness configuration consistently demonstrates a significant dampening effect on the downstream target tube across all flow velocities.Once the incoming flow velocity surpasses 0.2 m/s, there is a noticeable sharp increase in the amplitude of lower stiffness tube combinations, either upstream or downstream.This trend aligns with the characteristic behavior of fluid-elastic instability within tube bundles.The presence of lower stiffness tube combinations, whether upstream or downstream, to some extent, diminishes the critical velocity required for fluid-elastic instability, rendering the tube bundle more susceptible to instability.

Conclusions
The numerical calculation of coupling vibration of heat exchange tube bundles was completed using dynamic grid technology with rigid body motion equation.The vortex and vibration characteristics of tube bundles were analyzed, and the following conclusions were obtained: • Vortices shedding on the tube bundles are the primary cause of flow field fluctuation, which is the root cause of tube bundle vibration.The vortex shedding frequency on the bundles is the dominant frequency of the flow field.The tube bundles take the flow field dominant frequency as the primary vibration frequency and include the tube bundle's inherent frequency.• The vibration of multi-flexible tube bundles changes the vortex shedding distribution between tube bundles, and the vibration of the target tube differs from that of the single flexible tube bundle.On the premise that no resonance occurs, the upstream and downstream flexible tubes can suppress the vibration of the target tube to a certain extent, resulting in lower amplitude and more complex vibration frequency.However, the dominant frequency is the vortex shedding frequency.• The stiffness changes of the upstream and downstream tube bundles will affect the vibration of the central target tube.On the other hand, within the calculated range, no matter whether increasing the stiffness or decreasing the stiffness, the vibration of the central target tube cannot be better suppressed, and the amplitude of the central target tube is the smallest when the stiffness of the upstream and downstream and the target tube is the same.Reducing the stiffness of the tube bundles upstream or downstream will, to some extent, decrease the critical flow velocity of the target tube, making it more prone to experiencing fluid-elastic instability.This work can be used for better design and vibration suppression of shell-and-tube heat exchangers.

Figure 1 .
Figure 1.Profiles of failed tube bundles are shown in red rectangular (provided by Daqing Petrochemical Company).(a) Touching abrasion between tubes and baffle plate, (b) tube fatigue cracking.

Figure 2 .
Figure 2. Classification of heat exchange tubes.

Figure 3 .
Figure 3.The geometric model for turbulence model verification.

Figure 4 .
Figure 4. Zoom on the numerical mesh with scales.

Figure 5 .
Figure 5.Comparison between model calculation and experiment.

F 1 and
F 2 are mixed functions.According to Menter23 , the coefficients for the model are shown in Table2

Figure 6 .
Figure 6.Simplified vibration flexible tube model (m-mass of the tube; k-stiffness of spring; c-damping of the spring-damp-mass system).

Figure 7 .
Figure 7. Physical model of vibration of different stiffness tube bundles.

Figure 11 .
Figure 11.Time course curve of acceleration.

Figure 14 .
Figure 14.Amplitude time-frequency analysis of target tube under different stiffness combinations; (a) time domain; (b) frequency domain.

Figure 21 .
Figure 21.Amplitude time-frequency analysis of target tube under different stiffness combinations; (a) time domain; (b) frequency domain.

Table 2 .
Coefficients for the SST k-ω model.

Table 3 .
The mass, damping, and stiffness.-dimensional model.At the same time, only a dozen tubes in typical positions are selected for targeted research, and the calculation amount and time are reduced under the premise of accuracy.

Table 4 .
Stiffness and inherent frequency of tube.

Table 5 .
Calculation boundary and solution settings.

Table 6 .
Technical parameters of acceleration sensor.

Table 7 .
Amplitude of target tube under different stiffness combination schemes.